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ABSTRACT 


A study into the higher mode dynamic plastic response 
of beams is reported here. Within the scope of infinitesi- 
mal deflections, exact solutions are obtained for generic 
symmetric and antisymmetric modes, using rigid-plastic 
analysis. In order to account for the influence of finite 
deflections, an approximate procedure is used. Also a nu- 
merical procedure which satisfies the field equations at 
different time instants is developed. Calculations are 
performed for the elastic perfectly plastic case, using a 
finite-element formulation nlong with a finite-difference 
time integration. Results of the different solutions are 
compared with experiments. 
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1. INTRODUCTION 


Since the work of Lee and Symonds [1] and of Conroy | 2] 
who treated the problem of dynamic loading of rigid, perfec- 
tly plastic beams, one can identify three distint approaches 
that have commonly been used to solve the problems of dyna- 
mic loading of structures in the plastic range, namely: nu- 
merical treatments, generally including elastic and plastic 
behaviour, analytical solutions which usually neglect’ the 
elastic strains and the development of approximate procedu- 
res that predict only the overall characteristics of the dy- 
namic response. 

The numerical approach gives the possibility of analy- 
Sing the detailed response of complex structures, involving 
both elastic and plastic behaviour (material non-linearities) 
and geometrical non-linearities such as large deflections. 
Usually these procedures can incorporate also strain harde- 
ning and strain-rate effects in the material behaviour. 

The numerical schemes available are based in a_ special 
discretization of the structure and a timewise approximation 
of the response. Depending on the methods used for each, 
one can consider three basic types which use: (i) finite- 
-differences in space and time [3,4], (ii) finite-elements 


in space and finite-differences in time [5,6] and (iii) fi- 
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nite-element discretization both in space and time | 7,8]. 

These methods allow not only the representation of very 
complex structures but they also provide the time history of 
the transient response. They are appropriate for use in the 
final stages of the design of a structure, when there is 
enough detailed information to provide the necessary imput 
to these numerical procedures and, when it is required _ to 
verufy and adjust the details of the structural response. 

However, in the initial phases, since there is not 
enough information about the structure and since one is not, 
in general, interested in the detailed response, Simpler me- 
thods that provide a description of the overall behaviour, 
become more attractive. 

The main simplifications involved in the analytical — 
work consist in neglecting the elastic strains, by conside- 
ring the material to be rigid for stresses smaller than a 
given yield stress. This approximation is valid when the ki- 
netic energy of the disturbance is large compared with the 
energy that can be stored elastically in the structure. It 
has been shown [9,10] that the rigid perfectly plastic theo- 
ry 1s a reasonable first order theory as long as the ratio 
of the mentioned energies is at least greater than three. It 
was also observed [10] that for energy ratios greater than 


about ten, elastic vibrations do not have much effect on 


LZ 





the results. However, the duration of dynamic loading: should 
be short compared with the natural period of vibration of 
the structure. Since most of the problems of dynamic loading 
of structures in the plastic range are produced by blast or 
impact loading and are related to the desire of having the 
structure absorb as much energy by plastic work as possible, 
before fracture occurs, very often these requirements are 
met. 

This simple rigid plastic theory has been applied _ to 
beams [11,12], frames [13], plates [14,15] ‘and shells [16], 
and reasonable agreement with some experimental results 
{10,17,18] was obtained. 

However, these solutions use a small displacement as- 
sumption, in which the equilibrium equations are enforced in 
the undisturbed configuration of the structure. This assump- 
tion imposes a limit on the validity of the rigid-plastic 
idealization, since related to high levels of energy one has 
large deflections. The main developments of the rigid-plastic 
theory consisted in considering these large deflections or 
geometry changes [19-21], the incorporation of strain-rate 
sensitivity | 22-25] and strain-hardening | 26,27] in the ma- 
terial behaviour. When these effects were included a very 
good agreement with experiments were observed |[ 20,28,29]. 


The third line of thought, recognizes that even with 
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the rigid-plastic assumption, exact solutions are very diffi- 
cult to obtain, except for very simple situations, because of 
the discontinuities involved in the material behaviour. The- 
refore, to avoid the direct solution of the problem, approxi- 
mate procedures were developed in order to bound the final 
response or to approximate the deforming shape of the struc- 
ture. These are respectively the Bounding and the Wode Ap- 
proximation Techniques. 

The basis for the bounding techniques was laid down by 
Martin {| 30] who determined a lower bound for the response ti- 
me and an upper bound for the final deflection of impulsively 
loaded rigid-plastic structures undergoing infinitesimal de- 
flections. Moralles and Neville | 31] determined a lower bound 
for the final deflection and, Wierzbicki extended Martin's 
upper bound to the case of finite-deflections | 32]. Martin 
[33] determined bounds for the case of time-dependent plastic 
behaviour and recently Symonds and Chon | 34] developed an 
upper bound on deflection, valid for finite deflections and 
taking into account a time-dependent plastic behaviour. 

The mode approximation techniques, as presented by Mar- 
tin and Symonds [ 35 | are based on the uniqueness proof by 
Martin [36] which in turn is an extension of an earlier work 
by Tamuzh [37]. 


Mode solutions were defined as being composed of two 
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independent functions, one of space and other of times 


ee Cet, a xa) alee) 


where a isa scalar. The structure is therefore reduced 
to a one degree of freedom system since the deformation sha- 
pe remains constant while the magnitudes change with time.It 
was shown [35 | Licence Times tunetion, T(t), is linear, 
which implies that in a mode solution the acceleration is 
constant. The choice of the mode shape for a given problem is 
not arbitrary but instead results from the minimization of 
the initial value of an interrated quantity that measures 
the difference in kinetic energy between the actual struc- 
ture and the mode solution. This provides an upper’ bound 
on the error involved. This procedure, besides it's use- 
fulness in the estimation of final deformations, provides 
also a better insight for the main features of the de- 
formation process. 

The mode approach first developed for time indepen- 
dent material behavior was extended later to rate-sensiti- 
ve material [38] and to wider types of material behaviour 
[39]. 

Ho [40] extended the concept to "changing mode solu- 
tions“ which are solutions without the space-time separation. 


This intended to account for geometry changes and for moving 
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or stationary discontinuity interfaces. Convergence for dy- 

namically admissible solutions was proven. It was shown that 
the changing mode solution would converge to the stationary 

mode solutions of Martin and Symonds. 

While the changing modes would correspond to the last 
few phases of the actual motion of the structure, the statio- 
nary mode solution was the last phase. 

This extension allows a better accuracy in the predic- 
tions, especially for situations where earlier phases of mo- 
tion predominate but, on the other hand, it begins loosing 
the simplicity of analysis which is in fact one of the major 
attractions of the stationary modes. 

Lee [39] extended the approach, originally for impulsi- 
ve loading, to include dynamic loading problems. Also elas- 
tic, viscous, rigid perfectly plastic and rigid viscoplastic 
behaviour were considered and it waS given a wider meaning 
to mode solutions than just an approximation technique. 


Lee defined “instantaneous mode response" of the form: 


U, (x,t) = a(t) u;(x,t) (ieee) 


(x,t) = 8(t) O, (x,t) (1.20) 


and used the principle of virtual work to define the instan- 


taneous mode solution of specific problems. 
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He defines instantaneous elastic response as the kine- 
matically admissible field in which the potential energy is 
less than in all neighboring kinematically admissible fields 
which possess the same total displacement. When the external 
loads are zero, what corresponds to the case of free-vibra- 
tions, then the motion is in a natural mode of vibration, of 


the forms 


us (x,t) =o, (x) T(t) (1.3) 


as the classical normal modes of vibration. 

For the rigid-plastic case, the same approach was taken 
and an analogy was made with the upper bound theorem of Li- 
mit Analysis [41]. Lee concluded that the true collapse load 
of limit analysis is the primary instantaneous mode, as de- 
fined by him. 

Therefore, one could expect that the initial motion of 
a structure would follow the primary instantaneous mode res- 
ponse if the external load is arranged in such a way that it 
is on or slightly beyond the limit load. The implications of 
this fact will be discussed later in relation to Jones work 
[47,48]. 

Recently Symonds and Chon [42] extended the mode appro- 
Xlmation technique to accomodate finite deflections. 


When applying this technique to problems of finite de- 
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flections, because there are travelling hinges, a stationary 
mode solution is not applicable. Symonds and Chon conside- 
red then a sequence of instantaneous mode solutions at cho- 
osen instants of time. To go from one instant in time to 

the next one a numerical iterative procedure was used, so as 
to satisfy the field equations at that specific time instant. 

Basically their mode solution satisfies all the field 
equations but does not necessarily satisfy the stated ini- 
tial velocity condition. This corresponds to neglecting the 
initial phases of the motion or, the changing modes as de- 
fined by Ho. 

However, food results were obtained and it was obser- 
ved that the shape of the structure in sucessSive mode forms 
changed slowly. 

The authors noted that the numerical work in the appli- 
cation of the mode approximation, when finite-deflections 
are conSidered is much greater than in the infinitesimal ca- 
se. It is however substantially less than for the exact so- 
lution. Again, when extending the technique, one faces’ the 
problem of losing the simplicity and consequently it becomes 
less attractive. 

Recently, the bounding techniques which have been ba- 
sed mainly in energy considerations and the mode approxima- 


tion techniques which satisfy equilibrium equations after 
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having a given deformation mode choosen, seem to be mer- 
ging together. 

One of the difficulties in the mode technique consist 
in the determination of a kinematically admissible displa- 
cement or velocity field. 

Various procedures have been developed for a systema- 
tic determination of the approximate mode shape | 35,etc | 
and, recently it has been shown that extremum and variatio- 
nal principles can be used to determine mode shapes [43,44]. 

Independently, Jones extended the work of Sawcsuk [45, 
46] on finite deflections of plates loaded statically to 
the case of dynamic loading, developing an approximate pro- 
cedure for beams and arbitrarily shaped non-axisymmetric 
plates [47] which was later extended to shells [48]. Jones' 
method consists on the application of the principle of vir- 
tual velocities to a kinematically admissible collapse me- 
chanism. 

The main assumption is on the choice of the applicable 
collapse mechanism or displacement field. Jones assumed 
that the shape of the displacement field due to dynamic 
loads is the same as the one developed for the corresponding 
static collapse load. It was implicit on the assumption that 
the profile would be the same throughout the motion or, in 


other words, the hinges would remain stationary, although 
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plastic zones are allowed to be time-dependent. The final de- 
flection was then obtained by equating the rate of change of 
the initial enersy to the rate of dissipation of energzy by 
plastic work in the collapse mechanism. 

Wide comparisons with experimental results have been 
carried out [| 28,29,47-51] and good agreement was obtained. 

If one looks at this method in light of the so called 
mode approximation techniques, one realizes that it assumes 
a stationary mode throughout the motion. In a way, this can 
be considered a direct extension of the simple mode soluti- 
ons of Martin and Symonds [35] to the case of finite de- 
flections. It should be noted that Jones‘ method keeps’ the 
Simplicity of treatment while Symonds and Chon's extension 
[42] is expected to give better results in the cases of si- 
enificant changes of the hinge locations. 

The meaning of choosing the static collapse mechanism 
as the kinematically admissible displacement field can be 
established by refering to the work of Lee | 39]. This would 
correspond to choosing the primary instantaneous mode as 
being a stationary mode throughout the motion. 

Recently Jones and Wierzbicki [52] studied theoretical- 
ly and experimentally the higher modal dynamic plastic res- 
ponse of beams, associated with the determination of defor- 


mation mechanisms capable of absorbing more energy than the 
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first mode. They obtained an exact solution for the first 
three modes undergoing infinitesimal deflections and used 
Jones' method [47] to determine the influence of finite-de- 
flections on the first modal response. 

The present work consists of an extension of their 
analysis. 

Within the scope of infinitesimal deflections, exact 
solutions are obtained for generic symmetric and antisym- 
metric modes. In order to account for the influence of fi- 
nite deflections and to assess the relative merits of dif- 
ferent approximate procedures, Jones’ method [47] and Sy- 
monds and Chon's method [42] were applied to the first three 
modes. Extensions to generic symmetric and antisymmetric 
modes are incorporated. 

Finally a comparison is done with a more compreensive 
elastic-plastic numerical treatment of Wu and Witmer [6] 


and with experimental results. 
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2. INFINITESIMAL DEFLECTIONS 


2.1 Formulation of the problem 


The present formulation basically follows the one of 
Jones and Wierzbicki [52]. 

A fully clamped beam of lenght 2hL, subjected to diffe- 
rent forms of impulsive loading, is analysed herein. It is 
assumed that transverse displacements are small so that 
strains remain infinitesimal. 

Using the Bernoulli-Euler assumption of plane cross- 
-sections remaining plane, the displacement field of the 
beam may be approximated by the middle plane displacements. 
Moreover, if transverse shear deformations and rotary iner- 
tia are neglected, the equilibrium equations for a beam sub- 


jected to an impulsive transverse loading becomes 
Q' - pW = 0 (2.1la) 
M' +Q = 0 C2 rb)} 
where the positive sense associated with these quantities is 
shown in Figure 1, and the symbols are defined in Notation. 


Combining both equations, one obtains the equation of 


motions 


Br get 
a (2.2) 


— 
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The material is assumed rigid, perfectly plastic, neglecting 
therefore the effects of strain hardening and strain-rate 
sensitivity. The constitutive relations are then very sim- 
ple, as in Figure 2. For stresses smaller than the uniaxial 
yield stress, Jo» the strains are zero; at stress levels of 
Jo the strains become indefinitely large, of the same sign 
as the stresses; stresses higher than the yield stress are 


not allowed: 


€= 0 if lol <g 
O 


(2. 3a) 
Sign €=sign o i lo | = J. (2.3d) 
bal 0, (2.3¢) 


the corresponding, moment-curvature relations exibit the same 
behaviour [52]. 

To complete the formulation, boundary conditions are 
imposed as appropriate for each case. 

The procedure followed is to postulate a kinematically 
admissible displacement field, introduce it in the equilibrium 
equations, use the constitutive relations and appropriate 
boundary conditions to obtain the solution. Then verify whe- 
ther the solution is statically admissible. A solution that 
satisfies all these requirements is an "exact" solution. 

In all cases considered here the velocity field is of 


the form 
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w(x,t) => (x) W(t) (2.4) 
where @ {8s an x-dependent function describings the shape of 
the mode and W is a time dependent velocity amplitude as- 
sociated with the motion. 

Substituting (2.4) in the equations of motion (2.2.), 


results ins 


MY = -p Wo (2.5) 

It should be noticed now that W is independent of x 

and $¢ 1S independent of time. Therefore one can solve first 
the x-dependent part to obtain the shape and later the time- 


-dependent to obtain the motion. 
2.2 First, Second and Third Modes 


These three modes have been analysed by Jones and Wierz- 
bicki [52]. They used the time-independent acceleration fields 
shown in Figure 3. These acceleration fields were introduced 
in the equation of motion, (2.2) which was then integrated 
twice spacially to obtain the moment distribution. Use of the 
appropriate boundary conditions allowed the determination of 
the acceleration and location of the hinges. 


The acceleration for the first mode was found to be 





W= = 5 (2.6) 
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For the second mode it resulted 


7 M 

W = ene 32 Coa) 
where 

n= 2-V2 (2.8) 


and NL is the location of the hinge. 
It was shown that in the third mode the acceleration at 


the two hinges is the same and is given by (2.7), where now 
N= yo = 1 (258) 


The permanent transverse displacements were obtained inte- 
grating equations (2.6) and (2.7) in time. The results for 


the first, second and third modes are respectivelys 


W 

am. 2 dk (Qe) 

H 12 

bis 2 sA 

ait= (/9-1)° & (2.11) 

H 6 

iy 2d 

=e (2.12) 

H 2 

where 29 2 D 

ws L eV, 21 

eee = 5 (2.13) 
M, oH a Me 


is a dimensionless initial kinetic energy and Vy is the peak 
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initial velocity at the interior hinge locations. 
2.3 Symmetric Modes 


Symmetric modes are defined here as the ones 


the mode shape satisfies the following: 
¢(x) = $(2L-x) 


From an analysis of the first and third mode 
L521, it is apparent that only two types of shape 


9. exist in symmetric modes. 


One is applicable in the first region and is 
: ge 
Region ls >. = n , os x s1, 


the corresponding boundary conditions are: 


at x= 0 M = M, 

= N = 2 
at x 1 M Mo 
at x= on O20 


in which 


(27) 


results 


functions 


given by: 


(2.15) 


(2 6a) 
(22165) 


(2.16c) 


Substituting (2.15) in (2.5) and using (2.16), always 


results in one equation of motion. 


The other type of shape function is: 


Ce tas i iy 
Region i: 9%, = —= —4+—— . —icl ‘ 
Ww, 1, =", sre) 
1 1 lel 1 ile-l 
where the ee region has the yet and ce 
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ted at it's extreme points nN, and N,;_, , and has the as- 


sociated velocity amplitudes W. ; Wei » aS indicated in 


Figure 4, 


The applicable boundary conditions ares 


ac. «x See M= +t Mo G2 toa) 
al x =n, M = + MO Cate) 
at xen. Q = 0 (2.18c) 
at x meilis Q = 0 C2 lod) 


Use of these conditions in conjunction with equations 


(2.5) and (2.17) results in one equation of motion for the 


“un region and in the location of the Gar hinge. It 
should be noted that in the last region, one has: 
n, = L (2.19) 


i 
and therefore these are enough conditions to determine the 
location of all the hinges and to obtain one equation of mo- 
tion for each region. 

It should be noted that equations (2.15,2.16) reduce to 
the first mode expressions [52], when ",=L , and therefore 
we can obtain immediatly the solution for the first region 


from equation (2.6): 





W= -6 0 (2.20) 
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The solution for the zs region, is obtained by subs- 
tituting (2.127) in (2.4). 


Use of (2.18) leads to 


W. e¢ W. . = W (P) 


(eens) . 
1 irl u W (220) 


where W is Pivemupy (2.c0). 
Introducing (2.20) in (2.22) gives the general expres- 


th 


Sion for the location of the i hinge: 


1 ie (2.23) 
Now, when we are considering the last region of each 


mode shape, equations (2.19) and (2.23) result in 


fees 2 Ty (2, 24a) 


MSN oe PE Ip (2. 24b) 
Equations (2.24) lead to: 


N59 = L-2%27y (2325) 


sucessive substitution of (2.23) in (2.25) leads to the ge- 


neral forms 
l= 0.0 (2.26a) 
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-~ 1.3% Nn Pd ae 
eae = ]L-1°2 1 1 OP acame ae yt C2e7 GD) 


where nis the number of interior hinges per half span of 
the fen-1) °" mode. 
When i= n-l , (2.26b) becomes: 


— ae 
a = e227 } 


1+(n-1) v9 


Substitution of (2.27) in (2.20) and (2.21) gives the 
seneral expression for the acceleration, valid in all regi- 


ons s 


M ; 
W = - — ice (natn (2e 2) 
u L 
and (2.26b) becomes: 


n | 1-1" 2 

—t = itlisl? 2 ee le (25298 ) 
L 1+(n-1//2 

n =. 0 (2.29b) 


It should be noticed that for n=l and 2, equation 
(2.29) gives the previously obtained [52] hinge location for 
the first and third mode 


2.4 Antisymmetric Modes 


The antisymmetric mode shapes satisfy the following re- 


lations 
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(x) = - $(2L-x) (230) 


In these modes there exists three types of shape func- 
tions o.. Two of them are given by equations (2.15) and 
(2.17) of the symmetric modes and, the third type corres- 
ponds to the region just adjacent to the midspan of the 


beams 





n 
n2xsLb easy) 


the applicable boundary conditions ares 


at x= P M=+ M, (23248) 
at x= L 4 M= 0 C2.92b) 
at x =" = ; Q=0 (2.32c) 


where the second condition results from antisymmetric con- 
siderations. 
Introduction of (2.31) in (2.5) and satisfaction of 
(2.32) results ins 
n=, -%en 
Ase) ies Soe (2.33) 
Introducing equation (2.23) for i =n in equation 
(2.33) results ins 
n =p, - (4 
ney 7 E> (G+ 1)V2n, (2.34) 


SucesSive substitution of (2.23) in (2.34) results ins 
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Nn .=L- (g+ i)%2n G2 ) 


when i =n-] this expression will five the location of 
the first hinge after the supports, as a function of mode 


numbers 
n 


1 j 
is an (2.36) 
L  14%2 (n- S) 


Introducing this in equations (2.20) and (2.35) will 


result ins 


& 9 [1 + (n- 4) %2)? (2.37) 


and 


20 (2.38a) 


n. Ba 
c i “= al i= 1,2,...n (2.000) 
’ 1472 (n- 5) 


th mode. 


which are the general results applicable to the 2n 
For n= 1 equations (2.37) and (2.38) recover’ the 
results of the second mode, equations (2.7) and (2.8) as 


obtained in [52]. 


2.5 Equations of Motion for Symmetric and Antisymmetric 
Modes 


As resulted from the previous analysis of symmetric and 


antisymmetric modes, the equation of motion that governs all 
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the modes is the sames 





5 (2.39) 


where et is the distance of the first hinge from the clan- 
ped end of the beam, and assumes different values, aS given 
by equations (2.27) and (2.36). 

It should be noticed that for a fiven mode, it results 
from (2.3a) that the acceleration is constant. Indeed, Mar- 
tin and Symonds [35 ]showed that for stationary mode soluti- 
ons, like the present one, the velocity varies linearly 
with time, which is the same as having a constant accele- 
ration. 

The initial conditions applicable to the cases conside- 


red here are, at t= 0, 
W= 0 (2.40a) 
W= V (2.40b) 


what corresponds to applying an initial velocity of trian- 
gular shape and maximum value Te to an undeformed beams 
The shape of the initial velocity imparted to the beam is 
equal to the corresponding mode shape as determined from 


the previous analysis,equations (2.15,17,31) and Reference 


[52]. 


byz 





Interratins (2.39) twice with respect to time and ap- 
Mies the initial conditions (2.40), results in the deflec- 


tion amplitude: 





Nes Ve t 2en) 
a4 


The final detlection occurs when the velocity vanishes. 


This will happen ats 


oo O 
t, “eM (25?) 


and substituting in (2.41), we will obtain the final deflec- 


tion amplitude: 


Wo U ve né 
7 9 So (2.43a) 


0 


which can be rewritten as: 


2 
) (2.43b) 


“f 
H 


ale 

N 

“o— 

| = 
jot 


where A is given by (2.13). 
As the mode number increases, the corresponding value 


Of aT decreases and, consequently, the final deflections 


decrease also. 


The final deflected shape of the veam is fiven bys 


DE, 





Wi. n, a 
Tea ee CX) ( 2.) 


where % (x) is the appropriate mode shape. 
tor the first three modes (2.43) gives the results 


oztained by Jones and Wierzbicki | 52]. 
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3. APPROXIMATE PROCEDURE FOR FINITE DEFLECTIONS 


3.1 Description of the Method 


The approximate procedure used here is based on the the- 
oretical study of Jones [47] mentioned previously. It con- 
Sists of assuming a kinematically admissible collapse me- 
chanism and then,by using the principle of virtual veloci- 
ties, equate the energy imparted to the structure to the 
energy dissipated by plastic work at the hinges. 

The procedure was developed for a general non axis- 
symmetric, initially flat, rigid perfectly plastic plate and 
a specialised form was presented for a plate which deforms 
into a number of rigid regions separated by straight line 
hinges. 

The balance of energy for the particular case discussed 
above can be expressed by [47,52]: 


f  (P-pw) w dA = 
A 1 


EE oe ee 


f (Nw-M) 6, ae (3.1) 
1 £ 


where the area A on the left hand side of the equation is 
the total area of the plate which is subjected to the exter- 
nal pressure P, 6 


across the jon hinge and £ is the length of the hinge 


; is the relative angular rotation rate 


which, in the case of the beam is it's width. 
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The energy dissipated per unit width of the beam is fi- 


ven byt 


n ‘ 
D= (Nw-M) 6, (3.2) 


i=l 
while the left hand side of (3.1) represents the work done by 
the external and inertia forces. 

Equation (3.1) can be specialized for the case of a beam 
under impulsive loading by setting P=0 and performing one 


integration across the width of the beam, resulting in: 


mes 


Sf -pww dx = 


= (Nw-M) 6, (3.3) 
x 1 


Ji 


In the present study we will therefore apply equation 
(3.3) to the different modes. 

The imediate consequence of dealing with finite-deflec- 
tions for structures axially restrained consists on the de- 
velopment of membrane forces which increase in magnitude as 
the deflection increases, and provides a considerable reserve 
of strength. 

In passing, it should be noted that great care should 
be taken when modeling structures as axially restrained sin- 
ce Jones showed for the static case | 53] that remarkbly small 
in-plane displacements at the supports can change the respon- 


se from that of an axially restrained beam, with considerable 
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reserve strengthening beyond the limit lond for moderate In- 
ternal deflections, to that of a freely supported beam with no 
increase in strength beyond the limit load. 

Neglecting the effect of transverse shear forces | 54]as 
is appropriate at least for the lower modes [52], we must 
deal with the interaction curve that relates the moment and 
axial force. 

The Johanssen yield criteria, sometimes known as_ the 
Square yield curve, is known to represent well to behaviour 
of concrete structures | 55]. 

The Tresca and the Mises-Hencky yield criteria repre- 
sent better the behaviour of metal structures, and while the 
second is more accurate and more suited for numerical work 
[6], the Tresca yield condition being piecewise linear is 
in general choosen for analytical work. 

However, Johanssen's yield condition is still much sim- 
pler to use than Tresca and, is has been shown [| 25,47,56| 
that even for metal structures, a Square yield curve cir- 
cunscribing Tresca's curve could be considered to provide a 
lower bound to the maximum permanent transverse displacements, 
while an inscribing yield curve, as shown in Figure 8, (with 
Sides 0.618 times smaller) could be considered to give an 
upper bound. 


Therefore, in the present study both yield conditions 


Shi 





will be used and the results compared. 
In all cases, the assumed collapse mechanism has the 
same shape as the corresponding, infinitesimal] deflection re- 


sult, and because of symmetry only half-span is analysed. 
3.2 Johanssen Yield Criteria 


Using this criteria, the yield condition and normality 
requirements will be satisfied when one has at the hinges 


= = + 
N No and M zfs Moe 


It should be noted that for moderate deflections, in 
the absence of longitudinal accelarations, equilibrium re- 
quires the axial force to be constant throusshout the beam 


[57] and, consequently in the present case equal to No: 


In solving finite-deflection problems, very often a cons- 
tant axial force is assumed | 19,25,etc |]. 
It can be shown that the left hand side of equation 


(3.3) for half a beam is given by: 
5 = Ak 

f -uwwdx=-+ WW (3.4) 
3 


for all modes. 
Now, for each mode, the rotation rates should be deter- 
mined and introduced in (3.2) to determine the dissipated 


energy which, when equated to (3.4) will result in the 
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equations of motion. 

This procedure has been used by Jones and Wierzbicki, 
[52] who applied it to the first mode, They obtained for the 
maximum deflections 


W 
H 





m= Le (iea/3) 77-1) (3,5) 


Here we will apply it to the modes higher than the first. 
3.2.1 Second Mode 


The velocity field is given by equation (2.4) where the 
mode shape, shown in figure 7b, is the same as for infini- 


tesimal deflections [52] 1 


region 1: 2 CES ei (3. 6a) 
region 2 1 b5= =. ns xsl (3. 6b) 


The corresponding rotation rates are : 


; W 


av x= 0 3) Sy a: (3572) 
at x =n t= (ge 7D) 
n (L-n) 


where n iS given by equation (2.8), from the infinitesimal 
deflection analysis [52]. 


The normality requirements are satisfied when 
at x= 0 N= N : M= M sits) 


at x =n N= N  , M =-M (3. 8b) 


7 





If (3.7) and (3.8) are introduced in (3.3) we will 





obtains 
oe 2 2 
Wy + Y We + 3 = 6 359) 
where 
lig =) eo (3.10a) 
V 
O 
O 
Wy, = 7 (3.10c) 
ae 
and 
2 iy 
v= 4S. Tem ee 


The initial conditions are, at ty, = 01 
W, = 0 Co Ae 


Wy = 1 (3.12b) 


Solution of (3.9) and use of (3.12), in conjunction with 
the fact that the final deflection occurs when W, = 0 gives 


the duration of motion: 


¥ Y 


(Sigs) 


LO 





and the final deflections 


72 
Wat = 2 { (1+. 6 ¥2-8) mao ts Gan 
2 V2 3 


Results are shown in Figure 10. When A (6 V2-8)/3 <<1 
equation (3.14) reduces to the corresponding infinitesimal 


case result [52], equation (2.11) here. 
3.2.2 Third Mode 


As shown in Figure 9c, the velocity field is given by 


(2.4) where [52]: 


region ls d, = 2 » Os x sh (Sy ian 
Wy L-x x-N 

region 21 D5 a. a So = x (3.14b) 
Wo L-n Ln 


The rotation rates (3.7) remain valid in the present 


case and, in addition we have: 


Ow 
elie © E Ses an (3.15a) 


where now, from the infinitesimal deflection | 52], " is gi- 
ven by (2.9). 


Also, in addition to the normality requirements (3.8) 


we haves 


Ate =n N= ay a HM (3.15b) 


41 


the resultant equation of motion is given byt 


where (3.10) remains valid and now, 


(3.17) 


eg BE De: 
r 


4-342 


This equation predicts the final deflections 





War = 
which occurs at 


tef = = tan (3.19) 


¥ ¥ 


Results are shown in Figure 10. When A(80-57 J2)/(6-9 V2) <<1l 
equation (3.19) reduces to the corresponding infinitesimal 


case result [52], equation (2.12) here. 
3.2.3 Symmetric Modes 


The collapse mechanism has the same shape as given in 
section 2.3. The axial force is Ne throughout and normali- 
ty is satisfied having at the support M = M, and in conse- 


cutive hinges M = + Moe 


Equation (3.4) remains valid. To complete the formula- 


tion, we must now determine the dissipated energy, for a 


ho 





° ° h e e +? 
feneric i% hinge, as shown in Figure 4. 


! . Ge Eee ; 
The rotation rate across the 1 hinge, 18 @&lven bys 


Chega oe 4 ) 
6. = QW ee 3203 ) 


] 


The location of the re hinge in the symmetric modes, is 


riven by (2.29) and, if this is substituted in (3.20a) we 


obtains 
a = ew bh lGaeiy (3.20b) 
L 
The dissipation function corresponding to the ae hinge is 
then given by: 
V2 Mw. Ww 
D, = ——2— (144 #) (14+ 42(n-1)] (3.22) 
L H 


where use was made of equation (3.2). 
In the first region, which includes the hinge at. the 
Support and the first interior hinge, the dissipation func- 


tion is fiven bys 


. Wr. W 
D, = Mot L +(1+¥2)(1+4 wn, (34.2238) 
Semeusing (2.29), 


MW 
Dh ip 





{1+(1 +¥2)(1+4 5) }Li+(n-1) V2] (3222p) 


At the ae hinge, which is located at x=L, the dissipa- 
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tion funetion is «iven bys 


D _-_ =— Ds ee 3) 


einec, only half of the rotation is included becnuse only 
half of the beam is considered, 


Mhe total dissipation for the eo) ae mode iS siven 





byt 
D= D, + (n-3/2)D. (3, 24a) 
or ; 
MS My W 
He fre(aty ¥) [14248 (n-1)}} [1+ (n-1) V2) (3.24) 
H 
Equating (3.24) to (3.4) results in the equation of mo- 
tions 
‘9 2 ye Lr¥2 (n-1 
ety Wet SET =O ed 
where 


eee tn) ) Lit (n-1) V2 | (3.26) 
rN 


The duration of motion is given by 


tee = 2 tan”? (Lge gat (3.27) 


and the final deflection is: 
+V¥2(n- + Z 1/2 
ell me erey (it ee) -1) (3.28) 
| {it¥2 (n-1) | 


Wy 





For n=l and n=2 this result reduces to equations 
(3.5) and (3.18) of the first [52], and third mode respec- 


tively. 
3.2.4 Antisymmetric Modes 


As can be observed in Figure 5, in these modes only 
two types of dissipation function exist. One is valid for 
all but the region closer to the support and the other type 
is valid for the first region, containing the hinge at 


the support and one interior hinge. It 1s given by: 





MW | 
fee= —2— fit(iry 2) 2412969415 (n- 4 (3.29) 
1 L H 2 2 
if use is made of equation (2.36). 
The rotation rate across the on hinge is given by 


(3.20a). Using this equation and (2.36), one can derive the 
th 





applicable dissipation function for the i hinges 
2¥2 M 
O W 1 
Dp = ye (1th F) [1tv2 (n- $)] (3.30) 
The total dissipation for the an’? mode is given by: 
D=D, + (n-1) D; (Soles) 
or: ; 
Mo Ms W L ; 1 
D=— (1+(1+4 @)(1+2$2(n- 5)]} L1+Vv2(n- yi (30315) 
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Equating (3.31) to (3.4) results in the equation of mo- 


ie On 1 
a levees = (3.32) 

where 
= a [i+2 (n-1/2)] [1+2J2(n-1/2)] (Bessy, 


the final deflection is fiven by: 


_ _l+V¥2(n-1/2 m  _1422(n-1/2) 41/2 b 
Wee att (C1 3 ie B(n-t/2) >> 1} (3.3 


and it occurs at 


tee = tan+ [ 


1 _1+2(n-1/2 : 
f 


Y 246 (en-1/2 (3.35) 


When n= 1 these equations reduce to the second mode re- 


Pits, equations (3.9) to (3.13). 


3.3 Tresca Yield Criteria - Moderate Deflections 


In accordance with this criteria, the interaction curve 
between axial force and moment is parabolic, as shown in 


Figure 8 and can be expressed by [9]: 


2 
EM - el (3.36) 
O No 


The normality of the strain-rate vector to the yield 
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surface is Satisfied by the associate flow rules 


Mog NG 
a = ore [1 | AN (3.37a) 
N_ € ON 
O 
MoekK 
_ +s = A for |N/= N (7 37>) 
N € : 
O 


where € and are the generalized plastic strain rates, 
the extensiona] strain rate and the curvature rate respec- 
tively. Using this criteria, Jones [47] derived the dissipa- 
tion function for the interior hinge of a simply supported 


beam deforming in the first modes 


Be 
aL 


i 
WwW 


For the case of a clamped beam, it becomes: 


D. 


Oe. 
WwW P F 
: My (1+3 ye 0. (3.39) 


Also in a clamped beam,one has for the hinge at the supports 


D. 
1 


ae 
W 
H 
These dissipation functions are valid for N<N 9° 


When N=N., it can be seen from equation (3.36) that 
the moment must vanish. When that happens, these equations 
are no longer applicable. It can be shown that for a clamped 


beam that occurs when W, = l. 
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The present analysis are therefore valid for deflections 
smaller than the beam thickness which are called here modera- 
te deflections. 

In all the cases only half of the beam is analysed be- 


caum conditions of symmetry or antisymmetry apply. 
Bey3.l First Mode 


The velocity field is given by (2.4) where the shape is 
the same as used by Jones and Wierzbicki [52] for the infi- 
nitesimal case and when uSing the previous yield criteria. 


It is shown in Figure 8a and is given bys 


ee OsxslL (3.415) 


tl 


the corresponding rotation rates ares 


at x0 g=-+ (3.42a) 
at x=L a= = (3.42b) 


Using the appropriate rotation rates in conjunction 
with (3.39) and (3.40) and substituting all in (3.3), re- 


Sults in the equation of motions 


6M 
W + 2 


2 6 M, 
: (3.43) 








which agrees with Jones aga when his equations of motion 
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are specialized for the case of impulsive loading. 
This equation can be rewritten in the non-dimensional 


forms 


W, +8, W = B, (3.44) 


where (3.8) holds ands 


2 - & 
le (3.45a) 
m 226 
ay seo) 


the initial conditions to be satisfied are given by (3.12). 


Solving this equation by sucessive approximations, one 
obtains when using the second approximation and condition 


or l2): 


2 2 Pee 
B, t B, t 3B .4t B. t 
W, = t, + = —— -- - + (1+ —— = + Ser ) (3.46) 
and ; 
. Be ty 3B.t, 385 té 
Wy, = 148, ty = Spy (+ + So) (3.47) 


the duration of motion is obtained from equation (3.47) when 
Satisfying the condition Way = 0. Substitution of the smal- 
lest root of the resulting equation in (3.46) will give the 


final deflection amplitude or, the maximum deflection. 





3.3.2 Second Mode 


It can be shown that for this collapse mode, described 
in 3.2.1, the dissipation function applicable to the interior 
hinges is still given by (3.39). Using this equation and 


(3.40), the equation of motion obtained is still (3.44) whe- 





re nows3 
By = - 42 fuera) 
32-4 
Bo =-¢ é (3.48b) 


3 ¥2-4 


the solution for this mode is obtained from (3.46) and 


(3.47), using. the new values of B. 


bao.) Third Mode 


The same velocity field of section 3.2.2 applies here. 
The dissipation function for the hinge at the support and 
the hinge next to it 1s, as in the second mode, given by 
equations (3.39) and (3.40). For the central hinge, equation 
(3.38) is applicable. 


The equation of motion is given by (3.44) withs 


- & 642-5 
B= 5 (3.49a) 
+ 342-4 ° 
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- 2 (3,49b) 
32-4. 


yous 


the final deflection is obtained from (3.46) and (3.47). Re- 
sults from all these modes are shown in Figures 8-10 com- 
pared with the solutions obtained using the previous yield 


criteria. 
3.3.4 Symmetric Modes 


The total dissipation in these modes is obtained by 
considering all the three types of disspation function 


equations (3. 38-3.40). It is given bys 


= - : . Oe Or 
D Dy + D, + (n-2) D. + D, : n> ee 04) 


1 


for the (2n-1)"" mode. 
more tie first mode (n=l) since the central hinge is at 


the same time the first after the support, we haves 


(3. 50b) 


D. is the dissipation function at the support and 


is given bys 


D, = (1- 2 ) = [1tV¥2(n-1)] (3.51) 


where use was made of (3.40) and (2.27). 


Bat 





ror the ious hinge the rotation rate is given by (3.20b). 


Using this in conjunction with (3.38) and (2.27) results in 


2 Q 
p. = (1+4 45) * av2 (atv2(n-1)] (3.53) 
il & 
ial L 
The in Oe hinge occurs at mid-span and the applicable dissi- 


pation for all modes except the first 18S given by 
i = Sto 
=2d. (3.54) 


Therefore, for all but the first mode, combining (3.50) to 
(3.54), one obtains: 
MOW (w2 | 
oo A {> | 2+V2(8n-9) }+{ 2+2V2(n-1) 1} [i+V2(n-1)], na2 (3.55) 
L 


Equating (3.55) to (3.4) results in 


Wy +B) Wy 7 BS (3,56) 
where 
8, = ¢ L2+¥2(8n-9)] L1+V2(n-1)] (3. 57a) 
6 2 
B, = - > Litv2(n-1)] (3.57b) 


The solution to this equation is obtained from equations 
(3.46) and (3.47) using the present values of BB. When 


n= 3 equations (3.57) reduce to (3.49). 
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3.3.5 Antisymmetric Modes 


The dissipation function for the an vh 


mode 18 given by: 
D = Do + D, + (n-1) D; (oe50) 


where now 


N 


Dy = (te 15) FLU W2(n-2/2)] (3.59) 
Ps . 

D, = (143 =) . (1+?) | 1+$2(n-1/2) ] (3.60) 
fe 9 

D, = (1+4 “2 2 v2 [1+f2(n-1/2)] Couem) 


Combining these equations will obtains 


M 





W 2 
D=—2 {| 2+f2(8n-5)] oa + [ 242 ( 2n-1)] } L142 (n-1/2) J (3862) 
L H 


Substituting (3.62) in (3.3) and using (3.4), will re- 
sult in the equation of motion, equation (3.56), where now 
8, = ¢ L2tv2(8n-5)] L1+V2(n-1/2)) (3. 63a) 


B =-2 [ 2+f2(2n-1)] Li+V¥2(n-1/2) ] (35630) 


the solution of the equation of motion is given by (3.46) and 


(3.47). 


by 


a ba 





3.4 Tresca Yield Criteria - Large Deflections 


The results obtained in the previous section are valid 
for deflections smaller than one beam thickness. For sreater 
deflections, the beam exibits large deflections, as defined 
here, in which case it behaves like a plastic string. This 
has been reported elsewhere |e... 19,47]. 

We will analyse here only the first three modes since 
it is not likely to have modes higher than these enter this 


range of deflections. 
3.4.1 First Mode 


The dissipation function applicable in this case, as 


siven by Jones [47] iss 


ry W . 
) ; 7? 1 (3.64) 


- W 
[pe = at me 7 93 


1 


Equating this to (3.4), taking (3.42) into account, one 


obtains the equation of motion: 


W, + y° W, = 0 (3.65) 
where 
yo == (3.66) 


The initial conditions for this phase of motion are gi- 
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ven by equations (3.34) to (3.36), when WwW, = l. 

The duration of first phase is given by the smallest 
root of the equation resulting from substituting W, = 1 in 
equation (3.46). Substituting the resulting value of time 
in (3.47) will result in the initial velocity for this phase 
of motion, W,.. 


Then the initial conditions are at ty, = ty, 


W, = 1 (3.67a) 


een. (3.67b) 


% #7 
Solving equation (3.65) subjected to these initial con- 


dition results ins 


We: 
Wye = + siny(ty-ty,) + cosy(ty-t,,; ) (3.68) 


% 


the final deflection occurs when the velocity vanishes ats 


Sp ee 
tup = ty, + tan ~ ( a ) (3369) 


and the final deflection iss 

= 1/2 

Was 

alta) (3270) 
Y 


Wag 
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3.4.2 Second Mode 


Here, equations (3.65) and (3.68) to (3.70) of the 
previous analysis, are still valid. However, equation 


(3.66) becomes now 


y2 = 2 12 
dh 3¥2-4 (3.71) 


nN 





3.4.3 Third Mode 


The analysis remains the same as for the previous modes 


except that now the different rotation rates, lead to 


y= *< (5+392) (3.72) 
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4. APPROXIMATE NUMERTCAI PROCEDURE FOR FINTRIY Dielhiect LONS 


In the case of infinitesimal deflections, considered in 
Reference | 52] and in section 2 we were able to obtain an 
exact solution with a stationary mode shape because the 
initial velocity was arranged so as to have exactly the same 
shape as the mode solution and, because only bending moment 
was considered, consistent with the infinitesimal range of 
deflections. Only when these two conditions are satisfied 
Simultaneously can an exact solution with a stationary mode 
be obtained. 

For example, changing slightly the loading systems to 
beea Uniform initial velocity, instead of triangular, will 
result in having one phase of motion in which two hinges 
leave the supports and travel inwards. Only when the two 
hinges meet together will the mode shape considered in Refe- 
rence [52] or equation (3.41) prevail. 

Therefore a solution only with this shape as being 
stationary would only approximate the motion and therefore 
would not be exact. 

The inclusion of finite deflections in analysing cases 
where there are axial restrains, changes also the nature of 


the problem and stationary shapes can no longer be considered. 
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Therefore here we will consider velocity lelds of the 


forms 


w (x,t) = W(t) (x,t) (4,1) 


where the shape of the mode is allowed to change with time. 


The equilibrium equations for finite deflections are 


[57]: 


M - (N on asec (4,2a) 
oS eat (4. 2b) 
Nn (4. 2c) 


in the absence of transverse and longitudinal pressure loa- 


ding, From (4.2a) and (4.2c), it results: 


M -Nw +uw=0 (4,3) 


It is apparent from section 3 that the consideration of 
the Johanssen or of Tresca yield condition do not change the 
general concepts involved and that in general Johanssen's re- 
lation leads to less complicated treatments. Therefore, in 
this section, only the square interaction curve will be used 


to obtain the results. 


4,1 First Mode 


As now we will allow the mode to change with time, we 
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will not restrict the interior hinge to be located at cen- 
ter of the beam. Therefore, we will formulate the problem 
for the entire beam, of length 2L, while in the previous 
analysis only half span was considered. 

We will however be restricting the solution in that we 
will not account for situations such as two hinges spreading 
from the center of the beam, leaving a rigid region between 
or, any other mode form, except the one which has only one 
hinge in the interior of the beam. 


The mode shape is then fiven byt 


region ls os (cnc) = 07s x 3" (4,4a) 


region 21%, (x,t) = as ns x Ss eL (44d) 


where now 1 varies with time and its time derivatives must 
be considered. It should be noted that this mode shape is 
not a symmetric one and is therefore applicable to any tri- 
angular initial velocity condition, with it's peak at x=n. 

Equation (4.2c) together with the yield criteria impo- 
ses that N= No throughout. 


The applicable boundary conditions ares 


at x = 0 M = My (4.5a) 
= NN => 

at x M M, (4. 5b) 
at x = 2L M= M 
O 


De, 





at x =n Voss <0 (4, 5d) 
where V is defined as the vertical force and is given by: 
V = QtN) w (4.6) 
Introducing now (4.1) and (4.4a) in (4.3) and using 


conditions (4.5a,b,d) results ins 


7 oat 3 >N 6M 
i Se ih eco a = 5 = (0 (4,7) 








Use of (4.4b) and (4.5b-d) in the equilibrium equation, 
(4,3) results ins 


: Ww n a 6 Ws 





+ aes Eee aan (4,8) 
Diba) Pie) u(2L-" ) 


If one now introduces W from equation (4.8) in (4.7), 


the result iss 


2 6M 
n . (L)* +2 Ny (4.9) 
2u-" ywWwOL H 


this equation gives the velocity of propagation of the hinge. 
It should be noted that it vanishes when jn= L. 

The applied velocity field is symmetric and has it's ma- 
Ximum value at x=L. Then, the initial conditions are at 


t= O18 
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he oly (4,10a) 


W = 0 (4,10b) 


W = Vo (4,10c) 


Therefore the hinge is initially at "= L and the velo- 
city of it's propagation clearly vanishes. In this case we 
have therefore a stationary hinge as piven by this formula- 
THON 


Therefore equations (4.7) and (4.8) reduce just to 


. 3 No 6 Mo 
W + 5 TW + ae 0 (4,114) 


yp L uL 








It's non-dimensional form is simply: 


6 
A x = 0 Cola) ) 


It should be noted that because of the stationarity of 
the hinge, this solution reduces to the case considered by 
Jones and Wierzbicki [52] when taking finite deflections in- 
to account. 

Also, the stationarity of the hinge is the result of the 
symmetry of the initial velocity. The formulation, equations 
(4.8) and (4.9) handles the more general case of non-symme- 
tric initial conditions, in which case equation (4.9) pre- 


dicts a moving hinge. 
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4.2 Second Mode 


The velocity field is given in this case by (4.1) and 


(3.6) where in the last one " varies with time. The bounda- 


ry conditions are given by: 


at x= 
at x= 
at x= 
at x = 


0) M = 
n M = 
L M = 
n v= 


ikea) 
(4.12b) 


(cee Ze) 
(a2) 


where (4,12c) results from antisymmetry and V_ is given by 
(4.6). 

Introducing the velocity field corresponding to the 
first region in the equilibrium equations, (4.3) and using 
the appropriate boundary conditions results int 


3 N 6 M 


ms 2 =0 (4,13) 
N 


w—aWw- + 
n 








th 


- 


Performing the same for the second region, results ins 


a , 3.N 3 M 
Peco gy eae aay tt + 5 = 0 (4,14) 
L-1 (ttl ) u(L-") 


Now we must solve these two equations for the deflec- 
tion W and for the hinge location yn. In non-dimensional 


form they will be respectively: 
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n : 6 : ews : Mae — 
“on, 7 OR, * ‘el 
ee BC Zeis) az Wy 
ee (4,16) 
AN (1-7) X n,(1-",) 
where 
n n 
* => = (4,.17a) 
L 
1 ony 
e = v, 7 COE esos 


These equations will be solved numerically, as indicated la- 


ter. 
4.3 Third Mode 


The velocity field is given by (4.1) and (3.14) where 
N is to be understood as time varying. The applicable boun- 


dary conditions ares 


at x = 0 M = Mo (4,.18a) 
at x=" a (4.18b) 
at x=L ny eh (4,.18c) 
aa =o) Vo © (4,18d) 
at x=L V=0 (4.18e) 
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Introduction of the velocity Vtield (3.14) in the 
equilibrium equations, (4.3), and use of conditions (4b, U8n-c) 
leads to the expressions of moment in the two regions. Use of 
conditions (4.18d,e) will result in the equations of motion, 


which can be written as:3 








8 6 
Wo = oo LW, (1-3) -2& Wy) OAD) 
¥71 x nye Cae #oO ¥1 Noten 
W. +W oon n 
ve 8 #oO #7] eH ° @ He ae 
*0 ae @er ye fo aes ele (4,20) 
mM Won 12 W 
Ne — = + ° ~ r 2 (4.21) 
We Soom, We 0 Ue 


where now the velocities in the central region W, and on the 
lateral regions W are no longer equal. 


These equations will be solved numerically. 


4.4 Numerical Timewise Solution 


While for the infinitesimal case it was possible to 
obtain an analytical solution of the equations of motion, 
now we will solve them numerically since they don't seem 
amenable to a closed form solution. 


There exists a wide range of integration operators 
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that have been used sucessfully in structural dynamic pro- 
blems [58]. Their characteristics as to stability, converre 
and damping are well known [59], except for the recent and 
apparently more attractive Park's method [60,61]. 

Because in the present problem the elastic vibrations 
are neglected no major stability problems are expected and, 
instead of the more sofisticated operators such as Park and 
Houbolt , a Newmark ® operator [62]will be used. 

The Newmark £8 operator, for y= 1/2 and B = 1/4, also 
known as the trapezoidal method or the average acceleration 
method [63], exibits very good characteristics for linear 
problems, such as being unconditionally stable (B = 1/4) and 
introducing no false damping (y = 1/2). 

However, for non-linear problems, stability is not 
assured. 

The average acceleration method is an implicit method 
in that it needs information about the time step to be cal- 
culated, before it calculates it. Basically it uses for the 
value of acceleration during one time step the average value 
of the accelerations at its begining and end. Therefore, if 
no extrapolation scheme is used, aS was choosen here, itera- 
tion becomes necessary. 


In general, for the ee time step, the values of the 
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velocities and deflections are fiven bys 


W = 
n 


W 
n 


h 


where n 1s 


h 
e n ee oe 
W at oD (Wt) (4, 22a) 
h n e e 
Wiese ie (4,220) 
the nvh time step. 


For the present equations, one is dealing also with the 


hinge location. Considering the integration of this quantity 


and rewritting the equations (4.22) in a form more suitable 


to iteration, 


W 
n 


where 


the values of 


one hass 


role nlH 


A) To 


role 


nile 


Cy 23) 

(4,24) 

(4,25) 
h_ (4, 26a) 
ae (4.26b) 
oo (4, 26c) 


the acceleration Pe and hinge velocity ie 


are obtained from the equations of motion. 
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For the first mode, the hinge velocity was always zero 
and there is only one equation to be solved, (4.11), which 
turns out to have closed form solution [52]. However, it 
was also integrated numerically to verify the accuracy of 
the numerical calculations and, a perfect agreement was 
achieved. 

For the second mode we have to solve equations (4.15) 
and (4.16). The procedure used in this case was as follows: 


n 
to obtain the values at time step 1 = = hs begin by deter- 
t=] 


mining the constants (4.26) which will remain unchanged du- 


ring the iterations in the given time step; next assume that 


ey and determine (4.23) to (4.25); 


W and a = 


n-~ Wel 
substitute the calculated values sucessively in (4.16) and 


in (4.15) to obtain Ww, and "13 if these values differ from 


the initially assumed values by differences greater than the 
established convergence criteria, repeat the process until 
satisfactory convergence is attained; when the solution con- 
verged can go to next time step. The solution is considered 
to have converged when it agrees in the fourth decimal place. 
The procedure for the second mode is started by the initial 


conditions s 


W, = 0 (ERO 7) 
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Sel (4.27b) 


Oly 2- 2 (4.27c) 


The initial value of acceleration and hinge velocity is 
obtained by substituting (4.27) in (4.16) and (4.15). 
For the third mode three equations need to be solved, 


for W., W andn , (4.18-20). The general procedure is the 


1? 
Same as for the second mode except that now at each time 
step we need to achieve convergence in three variables 
instead of only two. 

In Figures 11 and 12 it is shown the time variation 
of all these variables for the first and second modes 
for a value of i» = 50. It should be noted that the equa- 
tions solved were in non-dimensional form and the results 
are applicable to the’ non-dimensional variables as defined 
before. 

Appendix contains a listing of the computer programs 
used to perform these computations. 

Calculations were made for different values of A so 
as to determine the variation of final deflection with A 
and these results are compared in Figures 8 and 9 with 


the previous solutions and with the experimental results. 
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The calculations made for the elastic-plastic behaviour 
of the beams, utilized the computer program JET 3, developed 
by Wu and Witmer [64]. 

The program, originally developed to calculate large 
elastic plastic dynamically induced deformations of free and 
restrained, partial and / or complete structural rings, can 
be used to analyse beams since it handles rings of arbitrary 
curvature. 

The program is based on the theoretical development re- 
ported in [6].It uses the assumed displacement version of 
the finite-element method to model the structure, together 
with the temporal finite-difference approximation. The equa- 
tions of motion for the spatial finite-element formulation 
are derived from the Principle of Virtual Work and D'Alem- 
bert’'s Principle. The Mises-Hencky yield criteria and it's 
associated flow rule are adopted to describe the elastic 
plastic behaviour of the material. The strain-hardening be- 
haviour is taken into account by using the mechanical sub- 
layer material model. The strain-rate effect is approxima- 
ted by assuming that the uniaxial stress-strain relation is 
affected the strain rate only by a quasi-steady increase in 


the yield stress above the static-test yield stress. The 
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program incorporates the Bernoulli-Euler ( or Kirchoff ) 
hypothesis, excluding transverse shear deformation. 
The behaviour of each finite-clement is characterized 


by four peneralized displacements at each node? 


where v and w are mid-plane displacements in the circun- 
ferential and normal direction, R the radius of curvature 
and the length cordinate. The displacement behaviour is 
represented by a cubic polynominal in n for both the circun- 
ferential, v , and the normal displacement w. 

Besides elastic restrains JET 3 includes three types 
of nodal displacement conditions namely, symmetry (v=y=0), 
ideally clamped ( v = w=yp = 0 ) and smoothly-hinged 
iy =w= 0 ). 

As forcing function, the program accepts initial veloci- 
ties and transient externally applied loads. These can _ be 
described as concentrated at the nodes, uniformly distribu- 
ted or with a Sine shape over the elements. 

The equations of motion are solved by applying an appro- 
priate timewise finite-difference operator which provides a 
solution step-by-step in finite-time increments. Two options 


are available for the integration operator: (1) the explicit 
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3-point central difference operator and Cy bie ting ek b 
Houbolt operator which uses a 4-point backward difference. 

The program is subdivided in four subprograms JET 3A 
to JET 3D which handle four different groups of problems. 

For the present calculations, it was choosen to use 
the Houbolt operator, Since it allows the use of larger ti- 
me steps and we are interested in the final deflection of 
the beams analysed. 

However, while for the central difference operator 
there exists a stability criterion for linear dynamic 
systems, which allow the estimation of the appropriate At 
to be used, no such criteria exists for the Houbolt opera- 
tor and the choice of At must be based on numerical ex- 
perimentation. 

Based on the numerical calculations reported in [6] 
for an impulsively loaded beam, it seemed appropriate to 
choose a St of Susec. 

While tor the central difference method unstable so- 
lution show unreasonably large deflections, for the Houbolt 
Operator a gradual degradation of the response may occur 
for large-deflection, non-linear response problems. There- 
fore, initially a calculation was performed also using A t= 
=3usec to assure thata At of 5usec was providing the 


converged solution. This was done for specimen 9 of ref. 
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| 52] which deforms in the second mode and identical results 
within plotting accuracy were obtained, aS Shown in Figure 
14. The rest of the calculations were performed using a At 
of 5psec. 

Because of antisymmetry considerations only half of 
the beam was analysed and 10 elements of equal length we- 
re used. The boundary conditions used were fully clamped at 
the support and smoothly hinyred at midspan. 

Basically the information the program needs consists 
of the mode locations and the shape of the prescribed ini- 
tial velocity. The initial velocity was given by it's value 
and slope at the nodal stations. As a result, the vicinity 
of the peak velocity was approximated by a higher order 
curve, resulting in higher velocities being introduced. The 
kinetic enersy that was given to the program resulted 2% 
higher that in the experiment. 

However, this would correspond to a higher velocity 
and a higher value of dA. In Figure 8 the final deflection 
is shown, corrected for new value of A, and Figure 14 
contains the full time history of the response. 

The final deflection is estimated by an average bet- 
ween it’s maximum and minimum values after some vibratory 


cycles have elapsed. 
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The same caleuluatbions were performed for specimen 10 
of Reference [52]. Now 1° elements of different size were 
used. A smaller one was centered at the peak velocity lo- 
cation and in it's vicinity finer mesh was used. The same 
type of result was obtained, as can be seen in Figures €& 


and 15. 
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6. EXPERIMENTAL DETAILS 


The general experimental arrangement used was similar 
to that employed by Jones in several previous studies | 28, 
29,49,50,52]. 

The present experiments represent a continuation of 
the studies reported in [52]. A reader is refered to it for 
amore complete discussion of the experimental details. 

The tests reported here were conducted on fully clamped 
beams made of aluminum 6061 16511 (E = 10.§x10° 1b/in’, 

P = 0.000251 lin acer The beams were nominally 5 in 
long 0.6 in wide and 0.2 in thick. They were loaded with 
layers of sheet explosive 3/8in wide, cut into diamond 
shapes, so as to produce a triangular impulsive velocity 
distribution. 

The impulse imparted to the beam depends not only on 
the weight of explosive but is also influenced by the leader 
arrangement and attennuator. Therefore calibration tests 
have been performed [52] to determine the necessary cons- 


tants. It resulted that for diamond shaped sheet explosive 


3/8in wide 


W, = W, + 0.03564 gm (6.1) 
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where Wo is the weirht of the explosive and, for the pre- 


gent explosive type the specific impulse iss 


J. = 0.2265 1b.sec/pm CGn2 ) 


The total impulse imparted to the beam is: 


Le ee (6.3) 
and the peak value of the velocity distribution is 


eee 
Vy 2 eee (6.4) 


where d is the length of the beam section which is covered 
with the exploSive. 

The ratio between kinetic energy and the maximum amount 
of strain energy which can be absorbed by a structure in a 


wholly elastic manner is estimated by: 


2 
0 EF ae 
B= (6.5) 
3 Oo 
+56 2 
where oO. = 49296 -1007 bya. 


The energy-absorbing characteristics can be measured 


by the dimensionless parameter 


W 
m 


= = seamen (6.6) 
Kv Mel. H~) 


WS 





where K) PomtnemcOtal initial kinetic energy. 


The permanent transverse displacements were given by 
the difference dial gauge readings before and after the test, 
with the beam still in the clamps. Only the maximum transver- 
se displacements, are shown in Table 1, and are plotted in 
Figures 8 , 9 and 1Q, together with results from Jones and 


Wierzbicki [52]. 
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7. DISCUSSION AND CONCLUSIONS 


A study into the higher mode dynamic plastic response 
of beams is reported here. Different solutions methods are 
used and results are compared with experiments. 

Tt is shown that the inelusion of finite-deflections 
Mmpreves Sipmificantliy the infinitesimal deflection solu- 
tions what agrees with previous results [19-21,25, pee 

The method developed by Jones (47] to account for 
the effect of finite-deflections proved simple to use and 
fave food results. When using the Tresca yield condition, 
as appropriate for metals, reasonable results were obtai- 
ned as can be seen in Figures 8 to 10. The use of a cruder 
description of material behaviour such as the square in- 
teraction curve, without strain-rate effects was much 
simpler to use and «ave once more | 47,51,52| surprizingly 
Boemuedesisn estimates of final deflection, equations 
em (3.14), (3.12%), as shown in Figures & to LO. 

The numerical procedure reported in section 4 which 
satisfies the field equations at every time step is more 
complicated and, for the present problem gave essentially 
the same results as using the previously mentioned method. 


it allowed for the hinge to propagate and it was found 
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bhat the hinge moved only slightly during the motion and 
the converrvence in the numerical iterative process was very 
ranid what confirms the observations of Symonds and Chon 
el? |. liowever, at each time step the moment was computed 
throuchout the beam and substantial yield violations were 
found to spread during the motion. 

The elastic-plastic numerical analysis is undoutly 
tne most compreensive and complex analysis method used. 
liowever, after being incorporated in a computer code, it 
becomes relatively simple to utilise. As can be seen in 
figure 9, 1t ives reasonable predictions, although with 
meremdency Of overestimating the deflections. This can 
also be observed in the results reported in [6] where the 
prediction of the final deflection of an impulsively loa- 
ded beam is somewhat higher than the experiment. 

It is apparent from the experimental results that, 
when the initial velocity is arranged with the shape of 
the mode considered, the beam deforms accordingly. Also, 
from experiments reported in | 52] with uniformly distri- 
buted load, it seem that the beams with different loading 
conditions tend to deform in the mode shapes reported. 
However, it 1s still to be known how a beam would deform 
when subjected to loadin; conditions substantially dif- 


ferent from the mode shape, Such as a triangular impulses 
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with it's peak not coinciding: with the maximum deflection 
of the mode shape. 

The use of normal modes 1S a powerfull tool when 
analysing the elastic vibrations of linear systems, the 
only case in which it is exact to use. towever, it has 
been used as an approximation in non-linear cases | 65-68]. 
1t rest to be determined if such concepts could be extended 
to plasticity and the higher dynamic plastic modes used in 


a modal superposition formulation. 
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APPENDIX 


Contains the results of the analysis described in 


section 4 together with the listing of the computer pro- 


grams useds 


Results for First Mode from Equation (4.11) 


Results for Second Mode from Equations (4.15,16) 


Program Listing: First Mode, Limited Interaction 


Program Listing: Seconde Mode, Limited Interaction 


OS 





WO A + 


<9 Gh ge 





1OA 





Suoz*TOT/ 


~taWVw 
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FIRST MODE+s LIMITED INTERACTION 


DIMENSION 8(200) 68D (200) *+BDD (200) +H(200) + T(200) »A4( 200) «8B (200) 


DIMENSION IT(200) 
DIMENSION F(2}) 
REAL MM(21) 


DEFINE CONSTANTS 
TWw=5 
RLAM=50,. 
DO 11 T=1,.200 
11 H(I)=0.01 
NREF=S 


DEFINE INITIat CONDITIONS 


0 
-6./RLAM 
PERFORM THE CALCULATIONS IN TIME STEPS 


WRITE (IWelfS) 
WRITE (IWe115) RLAM 


10S FORMAT (¢* s#teateeee FIRST MODE. LIMITED INTERACTION ###a4077) 


115 FORMAT (]X gst taeeaeee) AMDAHt,F7 Cy tettunvaranr sys //) 


1 CONTINUE 
N=Ne¢] 
NN=N=-1 

DEFINE THE CONSTANTS FOR FACH TIME STEP? 
TON) =T (NN) #H(N) 
A (NN) =BD (NN) *BND (NN) #H(N) Z20 
BB (NN) =B (NN) *BD (NN) #HOIN) 720 
IT(N) =0 
B(N) =R (NN) 

INITIATE ITERATION 
ITT=0 

S IT(NJ=IT(N) +1 
ITT=ITT+) 
00=B(N) 
BOD (N) =-6.%(1e42.*B(N) ) /RLAM 
BD (N) =A (NN) ¢+BDD(N) #H(N) S25 
B(N) =BB (NN) +BD(N) #H(N) 7/20 
IF (ITT eGT.e 5S) GO TO ]j 
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ERR=DD-B(N) 
IF (ABS(ERR) .GT. 0.0001) GO TO 5S 
M=N 
IF (N ~EQ. NREF) GO TO &80 
81 CONTINUE 
IF (N EQ. 200) GO TO SO 
IF (BO(N) LT. 0-05) GO TO 45 
GO TO ]l 
45 IF (BD(N) .LT. 0-000] ) GO TO SO 
H(N+#1)=0.0025 
GO TO } 


CALCULATE THE MOMENT AT DISCRETE POINTS ALONG THE BEAM 


80 NREF=NREF +s 
DEFINE THE POINTS TO BE CONSIDERED IN MOMENT CALCULATION 
FR=]./20- 
F(l)=0. 
00 2 1=292) 
Qo FCI) =F (I-11) +¢FR 
CALCULATE THE MOMENT 
00 40 J=1].2) 
MM (I) ==RLAM#BDDI(N) (F(T) #2 e-10) HFCL) S60 tl em le tF (T) 
40 CONTINUE 
WRITE (IWs1l10) N 
WRITE(IWs1N0) (F(T) s T=ls2)) 
| WRITE(IWs101) (MM(I) + T=ls2)) 
(100 FORMAT (1X.21F 6.2) 
(101 FORMAT (1X.21F6.2/) 
(110 FORMAT (1X, *##areNn=t, 3, tetaue) 
| GO TO a) 





MOTION FINISHED. PRINT THE RESULTS 
SO WRITE (IWs61) 
6] FORMAT (*)]?*) 
WRITE (IWs10S) 
WRITE (IWs115) RLAM 
i WRITE (IWs55) 
~~ 55 FORMAT (2X,* TIME 9X *DISPL og lIXs*VELPolIXs* ACCEL to SXotI Tt oSXsot*STE 
CP%e5X,*DELT®) 
| DO S6 N=lem 
— 56 WRITE (IWe60) T(N) sB(N) oBO(N) sBODI(N) o IT (CN) oNgHI(N) 
| 60 FORMAT (F639 3F 15 eS eSXeTOeSX oT 395K oF 50S) 
END 
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SECOND MODE. LIMITED INTERACTION 


DIMENSION &8(200) 6BD(200) +BDD(200) 6£ (200) +E D0 (200) +H (200) 
DIMENSION A(200) «BB (200) »C (200) e IT (200) sIST( 200) » T (200) 
DIMENSION F (50) oFF (50) 

REAL M1(50) -M2(50) 


DEFINE CONSTANTS 
IW=5 
RLAM=50. 
DO 11 I=1,200 
11 H(1)=0.005 
NREF=5 


DEFINE INITIAL CONDITIONS 


N=] 

T(1)=0. 

B(1)=0. 

BO(l)=l. 

E(1)=.5857864376 

BDD (1) =-3.8(22-E(N)) S(RELAMPE (1) FC LemMECL))) 
ED(1)=6./7 (QL AM#E (1)) *E (1) *BDD(1) 


CALCULATE TIME STEPS 


WRITE (IWel9S) 
WRITE (IWe115) RLAM 

1 CONTINUE 
N=Ne] 
NN=N-1 

DEFINE CONSTANTS FOR EACH TIME STEP 
T(N) =T CNN) #H(N) 
A (NN) =BD (NN) *BDD (NN) #H(N) 726 
BB (NN) =B (NN) *BD (NN) #H(N) 720 
C (NN) =E (NN) ED (NN) HON) 720 
IT (N) =0 
IST(N) =0 
B(N) =B(NN) 
E (N) =E (NN) 
BOD (N) =BDD (NN) 

INITIATE ITERATION 
ITT=0 

S IT(N) =ITI(N) «2 
ITT=ITT+} 
00=B (N) 
BO (N) =A (NN) *BDD(N) #H(N) Z20 
B(N) =BB (NN) BO (N) FH(N) 726 
IF (ITT eGT. 5S) GO TO } 
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aC BEGIN SUBITESATION 
| INN=0 
15 IST(N) =IST(N) +] 
INN=INNe 1] 
O=E (N) 
DB=BDN (N) 
ED(N) =6.%(1.*2.*%B(N) )/S(E (N) #BD(N) RL AM) +E (ON) BOD (N) SBD (N) 
E(N)=C (NN) +ED(N) HIN) 72, 
BDD(N) ==3.%(2e-E(N) +4.#R(N))S(QRLAM#E (N) #(1,-E(N))) 
IF (INN .GT.10) GO TO § 
ERROR=N-E (\)) 
IF (ABS(ERROR) .GT. 0.0091) GO TO 15 
ERRO=NB=-BDD (N) 
IF (ABS(ERRO) .GTe 0-000!) GO TO 15 
IC SUBITERATIONS FINISHED. GO TO NEXT ITERATION 
ERR=DD-B(N) 
IF (ABS(ERR) .GTe 0.0001) GO TO S 
M=N 
IF (N -EQ. NREF) GO TO R0 
81 CONTINUE 
IF (N .EQ2. 200) GO TO SO 
IF (BO(N) LT. 0-05) GO TO 45 
GO TO } 
4S IF (BO(N) .LTe 0.0001 ) GO TO 50 
H(N*1)=0.002 
| GO TO } 
| 80 NREF=NREF +5 
Ke CALCULATE THE MOMENT AT DISCRETE POINTS ALONG THE BFAM 


C DEFINE THE POINTS TO BE CONSIDERED 
FR=E(N) 7/20. 
F(1l)=0. 
DO 2 T=202) 
2 F(L)=F(I-1)+FR 
FF (1) =E(N) 
FFR=(1.°E(N)) 720. 
DO 3 T=2se)] 
3 FF CIT) SFR (I-1) +FFR 
C CALCULATE THE MOMENT 
DO 70 [=1sA1 
Mi (P)=(1.-208F (1) /E(N) ) “RL AM# (BOD (N) -BD(N) #ED(N) ZE(N)) #(F (1) ## 
C3-0-F (1) FE (N) F#2.00)/ (6,85 (N)) 
AM1=1.*E (N)-E(N) FE (N) 7/2. 
AMO=FF (I) *#3.-3e PF F(T) HC ote MFP (TL) HAM) +E CN) HHO, —2,%E (N) 
AM3=BOD (N) +BD(N) FED(N) / (Le sE(N)) 
M2 (LT) =—(leMFF (IT) SCL eo mE(N)) *RLAMPAM3BFAM2S (6.8% (1 2-E (N))) 
70 CONTINUE 
WRITE (IW9119) N 
WRITE (IW9100) (F(I)» T=lse2l) 
WRITE (IWe1lN0) (M1 (I). T=ls2)) 


a a 


C 


100 
101 
105 
110 
115 


WRITE (IWe 100) (FF (I) >. T=le2)1) 

WRITE (IWe101) (M2(T). I=le2l) 

FORMAT (1X.21F 6.2) 

FORMAT (1X.21F6.2/) 

FORMAT (* s#eeeteoee SECOND MODE. LIMITED INTERACTION ###at0 77) 
FORMAT (1X, ttteeen=l,[3Z, *Hetaste ) 

FORMAT (] Xo etttetatteet) AMDA=9,F 7 Ca ttttaettnnte ss ///) 

GO TO 81 


MOTION FINISHED. PRINT THE RESULTS 


90 
41 


55 


56 


60 


WRITE (IWs6)) 

FORMAT (1?) 

WRITE (IWs105) 

WRITE (IWe115) RLAM 

WRITE (IWe5S5) 

FORMAT (2X et TIME se OX eo tHINGE* glLOXes*DISPL* oll Xs*VEL toll Xe * ACCEL t's 9X, 


-*HINGEVEL # oSX ot LT te SXotSl* oSXeot*STEPt oOXe*DEL TT?) 


DO S6 N=lom 
WRITE (IWe60) TN) o£ (N) eBON) «BD (N) seBDDI(N) cED(N) oI TIN) sc ISTIN) oNeH(N 


C) 


FORMAT (F6,39SFI5S.Se7XeoTeOsSXeolT205XeoTl4eS5XeF 6.4) 
END 
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